%% report Field that must contain the following 
% 
% 3) *.stateSpectrum*  Spectrum 
% 
% 2) *.stateIRF* Impulse responses 
% 
% 1) *.stateMom* Will be usded for all moments (stds and cross-correlations)
% including variance decompositions, and also Historical Decompositions 
%  
% All of the above will be computed automatically for the observables 
% so no need to define them. 
% 
% =========================================================================

%% 1. Trim unfinished iterations 
% The unfinished iteration will have infs for their unfinished parts.  We
% can just grab the parts that had finished.
posTrim=isinf(abs(loopStru.logPost(:,2)));
if ~isempty(posTrim);
    loopStru.logPost=loopStru.logPost(posTrim==0,:);
    loopStru.logLikel=loopStru.logLikel(posTrim==0,:);
    loopStru.logPrior=loopStru.logPrior(posTrim==0,:);
    loopStru.parMode=loopStru.parMode(:,posTrim==0);
    loopStru.std=loopStru.std(:,posTrim==0);    
    loopStru.hessian=loopStru.hessian(:,:,posTrim==0); 
    loopStru.time=loopStru.time(posTrim==0); 
    loopStru.exitFlag=loopStru.exitFlag(posTrim==0); 
    loopStru.iterations=loopStru.iterations(posTrim==0); 
end
clear posTrim;
    
%% 2. *orderDescend* is the indicator of the ordered modes 
% Here, we grab the various statistics for the estimated iterations in
% order of their posteriors.
cd(cucd);
[~,orderDescend]=sort(loopStru.logPost(:, 2),'descend');

% Complete vector including calibrated parameters 
loopStru.parMode=loopStru.parMode(:, orderDescend);   
loopStru.logLikel=loopStru.logLikel(orderDescend,:);  
loopStru.logPost=loopStru.logPost(orderDescend,:);    
loopStru.hessian=loopStru.hessian(:,:,orderDescend);  
if loopStru.logPost(1,2) < max(loopStru.logPost(:,2) )
    error('Wrong Sorting of Results')
end
logPriorMat=loopStru.logPrior(orderDescend,:);
loopStru.std =loopStru.std(:,orderDescend);
optimMat=optimat(:,orderDescend);

%% 3. Assign model structure 
model.param=loopStru.parMode(:,1); 
model.savePath=location.savePath; 
model.handlePost=filterStru.funcPost; 
model.handleLikel=filterStru.funcLikel; 
model.logPost=loopStru.logPost(1,2);
model.logLikel=loopStru.logLikel(1,2); 

%% 4. Evaluate model at the mode 
%  Check posterior is correct 
fInitial=...
    feval(filterStru.funcPost,model.param(posStru.param.est),parStru,...
    posStru,model,dataStru,prpar,prss,filterStru,flags);
if abs(fInitial-loopStru.logPost(1,2)) >  10*estimopt.TolFun
    error('Log Posterior from funcPost and from minimization do not coincide')
end
clear fInitial; 

%% 5. Compute Hessian 
if flags.hessianDone==0
    if isfield(flags,'noHessian') && flags.noHessian
        stdHFD=nan( length(posStru.param.est) , 1 );
    else        
        disp('Begin Finite Difference Hessian');        
        [HFD,stdHFD,flagNotPD]=modeFinderHessian(estimopt.TolFun/10,...
            parStru,posStru,model,dataStru,prpar,prss,filterStru,flags);
        disp('End Finite Difference Hessian');
        crtcell([parnames(posStru.param.est) num2cprec([stdHFD(:) loopStru.std(:,1)]) ]);
        flags.hessianDone=1;
        close all;
        cd(location.savePath);
        save workspace;
        cd(cucd);
        disp('Hessian Done');
    end 
end

% ========================================================================
%% 7. Create ALL OUTPUT CELLS 
%% This is a mess that needs to be cleaned BIG time 
if ~flags.tableDone;
    modeFinderExcelTables;
    flags.tableDone=1;
end

%% 5. If ADD2FUNC exists, change function handle to this model after
if isfield(model,'add2Handle')==true
    if flags.initialCond==1
        disp('Cannot change handle if initial condition is used')
    else
        if isempty(model.add2Handle)==false && flags.handleChanged==0;
            model.handle=handleChange(model.add2Handle,model.handle,model.param,...
                model.solveOptions,model.addsol);
            flags.handleChanged=1;
        end
    end
end 
[GG, RR, CONS, eu, SDX, ZZ, initss, ssvec,~,names.steady,names.states,shockNamesStru]...
    =feval(model.handle,model.param,model.solveOptions,model.addsol);

model.CC=CONS;
model.RR=RR;
model.SDX = SDX;
model.ZZ = ZZ;

if isstruct(shockNamesStru)==true
    names.shocks=shockNamesStru.long;
    names.shockNames=shockNamesStru.short;    
else
    names.shocks=shockNamesStru;
    names.shockNames=names.shocks; 
end
[names.observables,posStru.obsInStates]=extractObsNames(ZZ(:,:,1),stateNames);
names.data=Ynames; 
names.param=parnames; 
clear GG RR CONS eu SDX ZZ initss ssvec; 
tempNames=fieldnames(prpar); 
for ii=1:length(tempNames); 
    fldnm=tempNames{ii}; 
    priorStru.(fldnm)=prpar.(fldnm); 
end 
clear tempNames; 

cd(location.savePath);
save workspace;


cd(cucd);
%% modelAnalysis
report.paramStd=stdHFD; 
if isfield(flags,'subSamples')== true && flags.subSamples==1 
    sampleStru.vector=sampleStru.sampleVec; 
    try
        sampleStru.cell=sampleStru.cell;
    catch
        sampleStru.cell=sampleStru.sampleCell;
    end
end 

% Version?
flags.OldMATLAB = verLessThan('matlab','8.4');

if flags.OldMATLAB
    modelAnalysisSummaryAJ_Old(model.handle,model.param,model,sampleStru,priorStru,report,location,names,posStru)
else
    modelAnalysisSummaryAJ(model.handle,model.param,model,sampleStru,priorStru,report,location,names,posStru)
end

%%
modeFinderStates; 